function varargout = myRdSRTM3(varargin)
%=========================================================================%
% Version : v1.1.1_20220619_alpha
% Date    : 20-Jun-2022
% 1. 新增一个功能
% method: 读取方式，'region'|'point'         
%         'method'不同前两个参数表示的意义不同：
%         'region'表示读取的地区坐标,此时显示LON和LAT之间的所划定的矩形范围
%         'point'表示点[lon1,lat1,alt1]到[Lon2,lat2,alt2]之间的所有高程。
% 语法：
% [a1,b1,c1] = myRdSRTM3([115.24,39.54,5],[116.24,40.34,5000],...
% 'ShowDEM','On','detailLevel',1,'method','point');
% [a,b,c] = myRdSRTM3([115.24,116.24],[39.54,40.34],...
% 'ShowDEM','On','detailLevel',1,'method','region');
%=========================================================================%
% Version : v1.0.2_20220607_alpha
% Date    : 9-Jun-2022
% 1.修改了一些路径问题
% 2.修复了输出参数的bug，之前的两个输出参数都是经度
%=========================================================================%
% Version : v1.0.1_20220607_beta
% Date    : 9-Jun-2022
% 1.注释和文档的问题
%=========================================================================%
% Version : v1.0.1_20220607_alpha
% Date    : 4-Jun-2022
% Author  : Xiaoyeyimier
% Email   : xiaoyeyimier@163.com
% Wechat  : xiaoyeyimier
%
% 读取大面积DEM文件（跨文件）
% 语法：
% (1) myRdSRTM3(LON,LAT);
% (2) myRdSRTM3(LON,LAT,'param','value');
% 'param':
% DEMPath: 数据库路径
% ShowDEM: 'On'显示地图;'Off'不显示地图
% detailLevel: 数值，表示压缩倍率。1时不压缩，小于1表示插值，插值使用三次插值。
%=========================================================================%
% 默认常数
GRIDSIZE = 1201; % SRMT3数据为1201×1201大小
GRIDWIDTH = 90;  % SRTM数据的采样间距为90m

%% 参数检查
if(nargin<2)
    error('myRdSRTM3 : at least 2 parameters required')
else
    % --可选参数 1数据库路径
    if any(strcmpi(varargin,'DEMPath')) % 参数对搜索
        tempIdx = find(strcmpi(varargin,'DEMPath'))+1;%名称-参数对中参数的索引
        if tempIdx>nargin
            error(['myRdSRTM3 syntax error : input variable ''DEMPath''',...
                ' must be formed through parameter-value pairs'])
        elseif ~isfolder(varargin{tempIdx})
            error('myRdSRTM3 syntax error : ''DEMPath'' must be a path');
        else
            CurrentDirection = varargin{tempIdx};
        end
    else
        CurrentDirection = cd;
    end
    % --可选参数 2绘图显示开关
    if any(strcmpi(varargin,'ShowDEM'))
        tempIdx = find(strcmpi(varargin,'ShowDEM'))+1;%名称-参数对中参数的索引
        if tempIdx>nargin
            error(['myRdSRTM3 syntax error : input variable ''ShowDEM''',...
                ' must be formed through parameter-value pairs'])
        elseif ~ischar(varargin{tempIdx})
            error('myRdSRTM3 syntax error : ''ShowDEM'' must be ''On'' / ''Off''');
        end
        switch lower(varargin{tempIdx})
            case 'on'
                ShowDEM = true;
            case 'off'
                ShowDEM = false;
            otherwise
                error('myRdSRTM3 syntax error : ''ShowDEM'' must be ''On'' / ''Off''');
        end
    else
        if nargout == 0
            ShowDEM = true; % 无参数输出、无特别指定画图时，自动画图
        else
            ShowDEM = false;% 有参数输出、无特别指定画图时，不自动画图
        end
    end
    % --可选参数 3缩放因子
    if any(strcmpi(varargin,'detailLevel'))
        tempIdx = find(strcmpi(varargin,'detailLevel'))+1; % 名称-参数对中参数的索引
        if tempIdx>nargin
            error(['myRdSRTM3 syntax error : input variable ''detailLevel''',...
                ' must be formed through parameter-value pairs'])
        elseif ~isnumeric(varargin{tempIdx})
            error('myRdSRTM3 syntax error : ''detailLevel'' must be a value');
        else
            detailLevel = varargin{tempIdx};
        end
    else
        detailLevel = 1; % 默认情况不压缩
    end
    if any(strcmpi(varargin,'method'))                % 读取方式的参数对
        tempIdx = find(strcmpi(varargin,'method'))+1; % 名称-参数对中参数的索引
        if tempIdx>nargin
            error(['myRdSRTM3 syntax error : input variable ''method''',...
                ' must be formed through parameter-value pairs'])
        elseif ~ischar(varargin{tempIdx})
            error(['myRdSRTM3 syntax error : input variable ''method''',...
                ' must be a string!'])
        elseif any(strcmpi(varargin{tempIdx},'region'))
            methodflag = 0;% methodflag 读取矩形区域
        elseif any(strcmpi(varargin{tempIdx},'point'))
            methodflag = 1;% methodflag 读取两点之间连线
        else               % 为其他模式保留
            methodflag = 0;% methodflag 读取矩形区域
            fprintf('myRdSRTM3 warning : ''method'' is forced to be region\n')
        end
    else
        methodflag = 0;  % methodflag 读取矩形区域
        detailLevel = 1; % 默认情况读取矩形区域
    end
    
end
% 更新本地数据库
DEMPath = fullfile(CurrentDirection,'DEMdat');% 下载文件存放在这里
MmapPath = fullfile(CurrentDirection,'m_map');% 下载文件存放在这里
addpath(DEMPath);% 添加路径
addpath(MmapPath);% 添加画图路径
batchProcessDEMFile(DEMPath); % 更新数据库
% 输入参数检查
if methodflag==0% 矩形模式
    [LON,LAT,scale] = ...
        LonLatCoordCheck(varargin{1},varargin{2},detailLevel,methodflag);
elseif methodflag==1% 点模式
    [Pt1LLA,Pt2LLA,scale] = ...
        LonLatCoordCheck(varargin{1},varargin{2},detailLevel,methodflag);
else% 为其他模式保留
end

%=========================================================================%
%% 多文件时，确定逐个读取文件的文件名
%% 矩形区域读取：
if methodflag==0
    % SRTM文件每个文件包含1°经度和1°纬度所覆盖的区域
    if max(LAT)==floor(max(LAT))
        lat = floor(min(LAT)):1:floor(max(LAT))-1;
    else
        lat = floor(min(LAT)):1:floor(max(LAT));
    end
    if max(LON)==floor(max(LON))
        lon = floor(min(LON)):1:floor(max(LON))-1;
    else
        lon = floor(min(LON)):1:floor(max(LON));
    end
    [matrixLon,matrixLat] = meshgrid(lon,lat);
    fileName = cell(length(lon)*length(lat),1);
    matrixSize = size(matrixLon);       % 还原数据时用[LATSIZE,LONSIZE]
    matrixLon = reshape(matrixLon,1,[]);% 变为一维数据方便不同模式处理
    matrixLat = reshape(matrixLat,1,[]);% 变为一维数据方便不同模式处理
    for fileIdx = 1:length(matrixLon)
        if matrixLat(fileIdx)<0
            stringLat = sprintf('S%02d',-matrixLat(fileIdx));
        else
            stringLat = sprintf('N%02d',matrixLat(fileIdx));
        end
        if matrixLon(fileIdx)<0
            stringLon = sprintf('W%03d',-matrixLon(fileIdx));
        else
            stringLon = sprintf('E%03d',matrixLon(fileIdx));
        end
        fileName{fileIdx} = [stringLat,stringLon,'.hgt'];
    end
    gsize = floor((GRIDSIZE-1)*scale);% 计算尺度放缩后的大小
    totalLat = zeros(1,gsize*matrixSize(1));
    totalLon = zeros(gsize*matrixSize(2),1);
    totalAlt = zeros(gsize*matrixSize(1),gsize*matrixSize(2));
    for DEMarrayIdx = 1:length(fileName)% 逐文件读取
        LonIdx = ceil(DEMarrayIdx/matrixSize(1));% 文件按照二维矩阵排布，计算其行索引
        LatIdx = mod(DEMarrayIdx,matrixSize(1)); % 文件按照二维矩阵排布，计算其列索引
        if(~LatIdx),LatIdx = matrixSize(1);end   % 列索引，第0列变为上一行最末列
        LonIdx = (LonIdx-1)*gsize+1 : LonIdx*gsize;% 数据点二维坐标，经度（行索引）
        LatIdx = (LatIdx-1)*gsize+1 : LatIdx*gsize;% 数据点二维坐标，纬度（列索引）
        % 先尝试读取本地文件，本地无对应文件时下载
        % 本地读取
        try
            tempDEMarrayStruct = readhgt(fullfile(DEMPath,fileName{DEMarrayIdx}));
        catch
            % 网络下载，现在应该不需要调整了（v1.0.1_20220607_beta）
            fprintf(['myRdSRTM3 WARNING : Cannot find DEM data %s in local database',...
                ',Trying to download it from the Internet...\n'],fileName{DEMarrayIdx})
            downloadLonRange = matrixLon(DEMarrayIdx);
            downloadLatRange = matrixLat(DEMarrayIdx);
            tempDEMarrayStruct = readhgt(downloadLatRange,downloadLonRange,'outdir',DEMPath);
        end
        % 放缩和文件合并
        if(DEMarrayIdx<=matrixSize(1))
            tempLat = tempDEMarrayStruct.lat(1:(GRIDSIZE-1));
            tempLat = resampleLL(tempLat,scale);
            totalLat((DEMarrayIdx-1)*gsize+1:DEMarrayIdx*gsize) = tempLat;
        end
        if(~mod(DEMarrayIdx,matrixSize(1)))
            tempLon = tempDEMarrayStruct.lon(1:(GRIDSIZE-1));
            tempLon = resampleLL(tempLon,scale);
            totalLon((DEMarrayIdx/matrixSize(1)-1)*gsize+1 : ...
                DEMarrayIdx/matrixSize(1)*gsize) = tempLon;
        end
        tempZ = tempDEMarrayStruct.z(1:(GRIDSIZE-1),1:(GRIDSIZE-1));
        tempZ(tempZ ==-32768) = 0;
        tempZ = resampleLL2D(tempZ,scale);
        totalAlt(LatIdx,LonIdx)= tempZ;
    end
    resLATIdx = totalLat>LAT(1)&totalLat<LAT(2);
    resLONIdx = totalLon>LON(1)&totalLon<LON(2);
    totalLat = totalLat(resLATIdx);
    totalLon = totalLon(resLONIdx);
    totalAlt = totalAlt(resLATIdx,resLONIdx);
%% 两点之间连线高程数据提取
elseif methodflag==1
    % 经纬高坐标检查
    Pt1ECEF = lla2ecef(Pt1LLA,'wgs84');
    Pt2ECEF = lla2ecef(Pt2LLA,'wgs84');
    LenPt1Pt2 = sqrt(sum((Pt2ECEF-Pt1ECEF).^2));
    directionVec = (Pt2ECEF-Pt1ECEF)/LenPt1Pt2;
    DirAcc = [(1:GRIDWIDTH/scale:LenPt1Pt2).'*...
        directionVec;(Pt2ECEF-Pt1ECEF)];% [(N*1)*(1*3);(1*3)]=>(N+1)*3
    PtNum = size(DirAcc,1);% 测试点数
    PtsECEF = repmat(Pt1ECEF,PtNum,1)+DirAcc;% 测试点坐标;
    PtsLLA  = ecef2lla(PtsECEF,'wgs84');
    PtsLL = floor(PtsLLA(:,[1,2]));% 获取所对应的SRTM文件（SRTM文件每个文件包含1°经度和1°纬度所覆盖的区域）
    fileName = cell(PtNum,1);% 文件名;
    for fileIdx = 1:length(PtsLL)% [lat,lon;lat,lon;...;lat,lon]
        if PtsLL(fileIdx,1)<0 % 纬度
            stringLat = sprintf('S%02d',-PtsLL(fileIdx,1));
        else
            stringLat = sprintf('N%02d',PtsLL(fileIdx,1));
        end
        if PtsLL(fileIdx,2)<0 % 经度
            stringLon = sprintf('W%03d',-PtsLL(fileIdx,2));
        else
            stringLon = sprintf('E%03d',PtsLL(fileIdx,2));
        end
        fileName{fileIdx} = [stringLat,stringLon,'.hgt'];
    end
    [fileName,PtIdx,~] = unique(fileName,'stable');
    PtsH = zeros(PtNum,1);
    for fileIdx = 1:length(fileName)% 逐文件读取
        if fileIdx~= length(fileName)
            tempPtIdx = PtIdx(fileIdx):(PtIdx(fileIdx+1)-1);
        else
            tempPtIdx = PtIdx(fileIdx):PtNum;
        end
        try
            tempDEMarrayStruct = readhgt(fullfile(DEMPath,fileName{fileIdx}));
        catch
            % 网络下载，现在应该不需要调整了（v1.0.1_20220607_alpha）
            fprintf(['myRdSRTM3 WARNING : Cannot find DEM data %s in local database',...
                ',Trying to download it from the Internet...\n'],fileName{fileIdx})
            downloadLonRange = PtsLL(PtIdx(fileIdx),2);
            downloadLatRange = PtsLL(PtIdx(fileIdx),1);
            tempDEMarrayStruct = readhgt(downloadLatRange,downloadLonRange,'outdir',DEMPath);
        end
        tempLat = tempDEMarrayStruct.lat;
        tempLon = tempDEMarrayStruct.lon;
        tempAlt = double(tempDEMarrayStruct.z);
        tempPtsLL = PtsLLA(tempPtIdx,[1,2]);
        PtsH(tempPtIdx) = interp2(tempLon,tempLat,tempAlt,tempPtsLL(:,2),tempPtsLL(:,1),'spline');
    end
    totalLat = PtsLLA(:,1);
    totalLon = PtsLLA(:,2);
    totalAlt = PtsH;
else
    %% 保留了其他读取模式的扩展
end
%% 输出参数控制
switch nargout
    case 0
        if ShowDEM
            MmapShowDEM(totalLat,totalLon,totalAlt);
        end
        return;
    case 1
        if ShowDEM
            MmapShowDEM(totalLat,totalLon,totalAlt);
        end
        dataOut.Lat = totalLat;
        dataOut.Lon = totalLon;
        dataOut.Alt = totalAlt;
        varargout = dataOut;
    case 2
        dataOut.Lat = totalLat;
        dataOut.Lon = totalLon;
        dataOut.Alt = totalAlt;
        varargout{1} = dataOut;
        if ShowDEM
            FigOut = MmapShowDEM(totalLat,totalLon,totalAlt);
        else
            FigOut = null(1);
        end
        varargout{2} = FigOut;
    case 3
        dataOut.Lat = totalLat;
        dataOut.Lon = totalLon;
        dataOut.Alt = totalAlt;
        varargout{1} = dataOut;
        if ShowDEM
            FigOut = MmapShowDEM(totalLat,totalLon,totalAlt);
            varargout{2} = FigOut;
        else
            varargout{2} = null(1);
        end
        Info.detailLevel = [num2str(detailLevel*100),'%'];
        dBytes = size(totalAlt,1)*size(totalAlt,2)*2;
        if dBytes>1e6
            dBytesStr = sprintf('%.2fMBytes',dBytes/1e6);
        elseif dBytes>1e3
            dBytesStr = sprintf('%.2fKBytes',dBytes/1e3);
        end
        Info.DataFootprint = dBytesStr;
        varargout{3} = Info;
end
end



%=========================================================================%
%% 私有函数，对一维的经度或者维度进行重新采样
function nDat = resampleLL(dat,scale)
if scale==1
    nDat = dat;
else
    delta = (dat(2)-dat(1))/scale;
    endOfVec = dat(end)+dat(2)-dat(1);
    if iscolumn(dat)
        dat0  = [dat;endOfVec];
    elseif isrow(dat)
        dat0  = [dat,endOfVec];
    end
    nDat  = dat0(1):delta:endOfVec;
    if abs(endOfVec-nDat(end))<delta
        nDat(end)=[];
    end
end
end
%=========================================================================%
%% 私有函数，对二维的高程数据进行三次插值
function nDat = resampleLL2D(dat,scale)
if scale==1
    nDat = dat;
else
    dat = double(dat);
    xGrid = 1:size(dat,1);
    yGrid = 1:size(dat,2);
    xNGrid = resampleLL(xGrid,scale);
    yNGrid = resampleLL(yGrid,scale);
    [XGRID,YGRID] = meshgrid(xGrid,yGrid);
    [XNGRID,YNGRID] = meshgrid(xNGrid,yNGrid);
    nDat = interp2(XGRID,YGRID,dat,XNGRID,YNGRID,'cubic');
    nDat = int16(floor(nDat));
end
end
%% 检查输入参数的完好
function  [Par1,Par2,scale] = LonLatCoordCheck(Par1,Par2,detailLevel,flag)
% 子函数用于检查数据的完好性
% 1.缩放倍率给出一定的提示
% 2.经纬度检查，经度在正负180度之间，纬度在正负90度之间。正表示东经/北纬
%   当使用矩形范围读取时，数据需满足Par1(1)<Par1(2);Par2(1)<Par(2)。
%   当使用两点连线读取时，数据需满足Par1和Par2都是1*3或者3*1经纬高坐标，且输
%   出会自动转为1*3行向量。
%-------------------------------------------------------------------------%
% 定义判断整数的匿名函数
isint = @(x) ~rem(x,1);
if detailLevel<1e-4                 % 控制压缩倍数以免过小
    detailLevel = 1e-4;             % 上限1/10000压缩
    fprintf(['myRdSRTM3 info : ',...% 
        'DEM data supports minimum 10000 times compression\n'])
elseif detailLevel>1e4              % 控制放大倍数以免过大
    detailLevel = 1e4;              % 上限10000倍放大
    fprintf(['myRdSRTM3 info : ',...% 
        'DEM data supports up to 10000 times interpolation\n'])
end
if detailLevel<1                    % 压缩时提示压缩倍率
    fprintf(['myRdSRTM3 info : ',...% 
        'Data is compressed by %d%%\n'],round(100-detailLevel*100))
end
scale = sqrt(detailLevel);          % 方法倍数
if ~isint(1/scale)&&~isint(scale)   % 插值导致原信息丢失时，提示
    fprintf(['myRdSRTM3 info : CUBIC INTERPOLATION was used in the',...
        ' output data, which may cause the original data lost!\n'])
end
% 检查画图范围
%-------------------------------------------------------------------------%
if flag==0
    if(length(Par1)~=2)||(length(Par2)~=2)
        error(['myRdSRTM3 input parameters'' error :',...
            ' Input LON/LAT must be a 1*2 array!'])
    end
    if Par1(2)<Par1(1)% 经度范围必须从小到大
        error(['myRdSRTM3 input parameters'' error :',...
            ' The first input longitude must be smaller than the second one']);
    end
    if Par2(2)<Par2(1)% 纬度范围必须从小到大
        error(['myRdSRTM3 input parameters'' error :',...
            ' The first input latitude must be smaller than the second one']);
    end
    if abs(Par1(1))>180||abs(Par1(2))>180 % 检查经度合理
        error(['myRdSRTM3 input parameters'' error :',...
            ' the absoulute value of the input longitude must be less than 180°'])
    end
    if abs(Par2(1))>90||abs(Par2(2))>90% 检查纬度合理
        error(['myRdSRTM3 input parameters'' error :',...
            ' the absoulute value of the input latitude must be less than 90°'])
    end
    if iscolumn(Par1),Par1=Par1.';end % 列向量转为行向量
    if iscolumn(Par2),Par2=Par2.';end % 列向量转为行向量
else
    if(length(Par1)~=3)||(length(Par2)~=3)% 检查数据大小=3
        error(['myRdSRTM3 input parameters'' error :',...
            ' Input Points must be in order of longitude, latitude and height'])
    end
    if abs(Par1(1))>180||abs(Par2(1))>180 % 检查经度合理
        error(['myRdSRTM3 input parameters'' error :',...
            ' the absoulute value of the input longitude must be less than 180°'])
    end
    if abs(Par1(2))>90||abs(Par2(2))>90% 检查纬度合理
        error(['myRdSRTM3 input parameters'' error :',...
            ' the absoulute value of the input latitude must be less than 90°'])
    end
    if iscolumn(Par1),Par1=Par1.';end% 列向量转为行向量
    if iscolumn(Par2),Par2=Par2.';end% 列向量转为行向量
    Par1 = Par1([2,1,3]);
    Par2 = Par2([2,1,3]);
end
end